!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Nonlinear, breaking mountain wave.
!!
!! \n
!!
!! This is the reference test case described in <a
!! href="http://journals.ametsoc.org/doi/pdf/10.1175/1520-0493%282000%29128%3C0901%3AAIOMPW%3E2.0.CO%3B2">[Doyle
!! et al., Monthly Weather Review, 2000]</a>.
!!
!! For the initial and boundary conditions, the experimental profiles
!! are used for \f$T\f$ and \f$u\f$, while the pressure profile is
!! computed by the hydrostatic balance:
!! \f{displaymath}{
!!  p(z) = p_s e^{\displaystyle -\frac{g}{R}\int_0^z
!!  \frac{d\zeta}{T(\zeta)}}.
!! \f}
!! Integrating the piecewise linear temperature profile we obtain:
!! \f{displaymath}{
!!  p(z) = p_k e^{\displaystyle -\frac{g}{R}
!!  \int_{z_k}^{z}
!!    \left[ T_k \frac{z_{k+1}-\zeta}{z_{k+1}-z_k}
!!         + T_{k+1} \frac{\zeta-z_k}{z_{k+1}-z_k} \right]^{-1}d\zeta
!!         } , \qquad z\in[z_k\,,\,z_{k+1}].
!! \f}
!! From this we have
!! \f{displaymath}{
!!  p(z) =
!!  \left\{\begin{array}{ll}
!!  \displaystyle
!!   p_k \left[
!!   \frac{(T_{k+1}-T_k)z+T_kz_{k+1}-T_{k+1}z_k}{T_k(z_{k+1}-z_k)}
!!   \right]^{\displaystyle 
!!     -\frac{g}{R}\frac{z_{k+1}-z_k}{T_{k+1}-T_k}}, &
!!   T_{k+1}\neq T_k \\[4mm]
!!   p_k e^{\displaystyle -\frac{g}{RT_k}(z-z_k)} \left[
!!   \displaystyle 1 +
!!   \frac{g(z-z_k)^2}{2RT_k^2(z_{k+1}-z_k)}(T_{k+1}-T_k) \right] , &
!!   T_{k+1}\to T_k.
!!  \end{array}\right.
!! \f}
!<----------------------------------------------------------------------
module mod_boulder_windstorm_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

 use mod_physical_constants, only: &
   mod_physical_constants_initialized, &
   t_phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_boulder_windstorm_test_constructor, &
   mod_boulder_windstorm_test_destructor,  &
   mod_boulder_windstorm_test_initialized, &
   test_name,        &
   test_description, &
   phc,              &
   ntrcs,            &
   ref_prof, t_ref,  &
   coeff_dir,        &
   coeff_neu,        &
   coeff_norm,       &
   coeff_sponge,     &
   coeff_visc,       &
   coeff_init,       &
   coeff_outref

 private

!-----------------------------------------------------------------------

 ! public interface
 character(len=*), parameter   :: test_name = "boulder-windstorm"
 character(len=100), protected :: test_description(2)
 type(t_phc), save, protected  :: phc
 integer, parameter            :: ntrcs = 1
 character(len=*), parameter   :: ref_prof = "isothermal"
 real(wp), target              :: t_ref

 logical, protected :: mod_boulder_windstorm_test_initialized = .false.
 ! private members
 real(wp), parameter :: &
   zu0(24,2) = reshape( (/ &
    ! velocity data
      0.00e3_wp ,  9.0_wp , &
      1.25e3_wp ,  9.0_wp , &
      2.50e3_wp , 16.0_wp , &
      3.75e3_wp , 30.0_wp , &
      5.00e3_wp , 35.0_wp , &
      6.25e3_wp , 42.5_wp , &
      7.50e3_wp , 38.0_wp , &
      8.75e3_wp , 32.0_wp , &
     10.00e3_wp , 31.0_wp , &
     10.75e3_wp , 38.0_wp , &
     11.25e3_wp , 30.0_wp , &
     12.00e3_wp , 34.0_wp , &
     12.50e3_wp , 17.0_wp , &
     13.00e3_wp , 17.5_wp , &
     14.00e3_wp , 30.0_wp , &
     15.00e3_wp , 28.0_wp , &
     16.25e3_wp , 21.0_wp , &
     17.50e3_wp , 15.0_wp , &
     18.75e3_wp , 13.0_wp , &
     20.00e3_wp ,  5.0_wp , &
     21.00e3_wp ,  0.0_wp , &
     22.00e3_wp ,  3.0_wp , &
     23.75e3_wp ,  1.5_wp , &
     25.00e3_wp ,  7.5_wp   &
    ! end velocity data
       /) , shape(zu0) , order=(/2,1/) ), &
   zt0(19,2) = reshape( (/ &
    ! temperature data
      0.00e3_wp , 272.0_wp , &
      0.50e3_wp , 273.0_wp , &
      2.50e3_wp , 256.0_wp , &
      4.00e3_wp , 250.0_wp , &
      5.00e3_wp , 246.0_wp , &
     10.00e3_wp , 208.0_wp , &
     11.20e3_wp , 209.0_wp , &
     11.25e3_wp , 214.0_wp , &
     12.50e3_wp , 216.0_wp , &
     14.00e3_wp , 208.0_wp , &
     16.25e3_wp , 213.0_wp , &
     17.00e3_wp , 220.0_wp , &
     17.50e3_wp , 220.0_wp , &
     18.75e3_wp , 217.0_wp , &
     20.00e3_wp , 221.0_wp , &
     21.25e3_wp , 223.0_wp , &
     22.50e3_wp , 222.0_wp , &
     23.75e3_wp , 222.0_wp , &
     25.00e3_wp , 220.0_wp   &
    ! end temperature data
       /) , shape(zt0) , order=(/2,1/) )
 real(wp) :: p0(size(zt0,1))

 real(wp) :: &
   nu,     & ! kinematic viscosity
   ! mountain geometry
   hm, a, x0, &
   ! sponge layers
   zd,    & ! begin of the upper damping layer
   zt,    & ! domain top
   alpha, & ! gravity wave damping
   dv,    & ! vertical thickness acoustic damping
   dh,    & ! horizontal thickness acoustic damping
   x1(2), x2(2) ! lateral domain boundaries
 character(len=*), parameter :: &
   test_input_file_name = 'boulder_windstorm.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_boulder_windstorm_test'

 ! Input namelist
 namelist /input/ &
   nu, hm, a, x0, zd, zt, alpha, dv, dh, x1, x2

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_boulder_windstorm_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, i, ierr
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)   .or. &
       (mod_kinds_initialized.eqv..false.)      .or. &
       (mod_fu_manager_initialized.eqv..false.) .or. &
       (mod_physical_constants_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_boulder_windstorm_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Windstorm test, dimension d = ',d
   write(test_description(2),'(a,e9.3,a)') &
     '  viscosity nu = ',nu,'.'

   ! All the physical constants use the default values: phc unchanged

   t_ref = 300.0_wp

   ! Initialize the module array
   p0(1) = phc%p_s
   do i=2,size(p0)
     p0(i:i) = pint(zt0(i:i,1),zt0(:,1),p0,zt0(:,2))
   enddo

   mod_boulder_windstorm_test_initialized = .true.
 end subroutine mod_boulder_windstorm_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_boulder_windstorm_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_boulder_windstorm_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_boulder_windstorm_test_initialized = .false.
 end subroutine mod_boulder_windstorm_test_destructor

!-----------------------------------------------------------------------

  !> Viscosity
  pure function coeff_visc(x,uu) result(mu)
   real(wp), intent(in) :: x(:,:), uu(:,:)
   real(wp) :: mu(2,size(x,2))

    mu = nu
  end function coeff_visc

!-----------------------------------------------------------------------

 !> Hydrostatic atmosphere
 !!
 pure function coeff_init(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))

  integer :: l
  real(wp) :: r2
  real(wp), dimension(size(x,2)) :: z, p, t, uuh

   z = x(size(x,1),:) ! vertical coordinate
   ! primitive variables
   t   = linint(z,zt0(:,1),zt0(:,2))
   p   = pint(z,zt0(:,1),p0,zt0(:,2))
   uuh = linint(z,zu0(:,1),zu0(:,2)) ! horizontal velocity

   ! translate the result in conservative variables
   u0(1,:) = p/(phc%rgas*t)                                     ! rho
   u0(2,:) = u0(1,:)*(phc%cv*t + 0.5_wp*uuh**2 + phc%gravity*z) ! e
   u0(3,:) = u0(1,:)*uuh                                        ! Ux
   u0(4:2+size(x,1),:) = 0.0_wp                                 ! Uyz
   ! subtract the reference state
   u0(1,:) = u0(1,:) - u_r(1,:) ! rho
   u0(2,:) = u0(2,:) - u_r(2,:) ! e

   ! tracers
   if(ntrcs.gt.0) then
     do l=1,size(x,2)
       r2 = (x(1,l)-x0)**2 + 0.25_wp*(z(l)-zd/2.0_wp)**2
       u0(2+size(x,1)+1,l) = 1.0_wp*exp( -r2/(zd/3.0_wp)**2  )
     enddo
   endif

 end function coeff_init
 
!-----------------------------------------------------------------------

 !> Use the initial condition.
 pure function coeff_outref(x,u_r) result(u0)
  real(wp), intent(in) :: x(:,:), u_r(:,:)
  real(wp) :: u0(2+size(x,1)+ntrcs,size(x,2))

   u0 = coeff_init(x,u_r)
   u0(2+size(x,1)+1,:) = 0.0_wp

 end function coeff_outref

!-----------------------------------------------------------------------

 !> Use the initial condition.
 pure function coeff_dir(x,u_r,breg) result(ub)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: ub(2+size(x,1)+ntrcs,size(x,2))

   ub = coeff_init(x,u_r)

 end function coeff_dir

!-----------------------------------------------------------------------

 pure function coeff_neu(x,u_r,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  real(wp), intent(in), optional :: u_r(:,:)
  integer, intent(in), optional :: breg
  real(wp) :: h(size(x,1),1+size(x,1)+ntrcs,size(x,2))

   h = 0.0_wp

 end function coeff_neu

!-----------------------------------------------------------------------

 !> Normal to the boundary
 !!
 !! \warning This function should be only called for the bottom
 !! boundary.
 !!
 !! In 2D, we have
 !! \f{equation}{
 !!  y(x) = \frac{h_m}{1+\left(\frac{x-x_0}{a}\right)^2}
 !! \f}
 !! so that
 !! \f{equation}{
 !!  \tilde{\underline{n}}(x) = 
 !!   \left[\begin{array}{c} \frac{dy}{dx}(x) \\\ -1 \end{array} \right]
 !!   = \left[\begin{array}{c} 
 !!    -2\frac{h_m}{a}\frac{\xi}{\left(1+\xi^2\right)^2} \\\ -1
 !!   \end{array} \right]
 !! \f}
 !! where the \f$\tilde{}\f$ indicates that the resulting vector must be
 !! normalized and
 !! \f{equation}{
 !!  \xi = \frac{x-x_0}{a}.
 !! \f}
 !<
 pure function coeff_norm(x,breg) result(nb)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: nb(size(x,1),size(x,2))

  integer :: l
  real(wp) :: xi(size(x,2))

   breg_case: select case(breg)
    case(3) ! bottom boundry, 2D
     xi = (x(1,:) - x0)/a
     nb(1,:) = -2.0_wp*hm/a * xi/(1.0_wp+xi**2)**2
     nb(2,:) = -1.0
     do l=1,size(nb,2)
       nb(:,l) = nb(:,l)/sqrt(nb(1,l)**2+nb(2,l)**2)
     enddo
    case default
     nb = huge(1.0_wp) ! this is an error
   end select breg_case

 end function coeff_norm

!-----------------------------------------------------------------------

 pure subroutine coeff_sponge(tmask,tau,taup,x)
  real(wp), intent(in) :: x(:,:)
  logical, intent(out) :: tmask
  real(wp), intent(out) :: tau(:), taup(:)

  real(wp), parameter :: pi_trig = 3.141592653589793_wp, &
    big = huge(1.0_wp)
  integer :: l, i
  real(wp) :: z, zn, d_bnd, coeff

   ! initialization
   tau = 0.0_wp
   taup = 0.0_wp
   tmask = .false.

   do l=1,size(x,2)

     ! upper layer
     z = x(size(x,1),l)
     if(z.gt.zd) then ! upper layer nonzero

       tmask = .true.

       ! gravity wave damping
       zn = (z-zd)/(zt-zd) ! normalized vertical coordinate
       if(zn.le.0.5_wp) then
         tau(l) = tau(l) + alpha/2.0_wp*(1.0_wp-cos(zn*pi_trig))
       else
         tau(l) = tau(l) + alpha/2.0_wp*(1.0_wp+(zn-0.5_wp)*pi_trig)
       endif

       ! acoustic wave damping
       d_bnd = (zt-z)/dv ! normalization
       if(d_bnd.le.3.0_wp) then
         ! tmask already set to .true.
         if(d_bnd.lt.1000.0_wp/big) then
           coeff = big/1000.0_wp
         else
           coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
         endif
         tau(l)  = tau(l)  + coeff
         taup(l) = taup(l) + coeff
       endif

     endif

     ! lateral sponges (only acoustic wave damping)
     d_bnd = big
     do i=1,size(x,1)-1 ! first d-1 coordinates
       d_bnd = min( d_bnd , min( abs(x(i,l)-x1(i)) , abs(x(i,l)-x2(i)) ) )
     enddo
     d_bnd = d_bnd/dh ! normalization
     if(d_bnd.le.3.0_wp) then
       tmask = .true.
       if(d_bnd.lt.1000.0_wp/big) then
         coeff = big/1000.0_wp
       else
         coeff = (1.0_wp-tanh(d_bnd))/tanh(d_bnd)
       endif
       tau(l)  = tau(l)  + coeff
       taup(l) = taup(l) + coeff
     endif

   enddo
   
 end subroutine coeff_sponge

!-----------------------------------------------------------------------
 
 !> 1D linear interpolation
 !!
 !! It is assumed that \c xn is sorted in increasing order. The first
 !! and last values are extrapolated as constants.
 !<
 pure function linint(x,xn,yn) result(y)
  real(wp), intent(in) :: x(:), xn(:), yn(:)
  real(wp) :: y(size(x))

  integer :: nn, i, k
  real(wp) :: fdx(size(xn)-1), q
 
   nn = size(xn)
   fdx = 1.0_wp / ( xn(2:nn) - xn(1:nn-1) )

   do i=1,size(x)
     if(x(i).le.xn(1)) then
       y(i) = yn(1)
     elseif(x(i).ge.xn(nn)) then
       y(i) = yn(nn)
     else
       do k=2,size(xn)
         if( x(i).le.xn(k) ) then
           q = fdx(k-1)*(x(i)-xn(k-1))
           y(i) = yn(k-1)*(1.0_wp-q) + yn(k)*q
           exit
         endif
       enddo
     endif
   enddo
 
 end function linint
 
!-----------------------------------------------------------------------

 !> 1D pressure interpolation
 !!
 !! This function is similar to \c linint, but uses the nonlinear
 !! expression of the pressure.
 !!
 !! \note This function can be used to initialize the array of the
 !! nodal pressure, provided that the first element is set (surface
 !! pressure) and that it is called for increasing nodal vertical
 !! coordinate.
 !<
 pure function pint(z,zn,pn,tn) result(p)
  real(wp), intent(in) :: z(:), zn(:), pn(:), tn(:)
  real(wp) :: p(size(z))

  real(wp), parameter :: dt_toll = 1.0e4_wp*epsilon(1.0_wp)
  integer :: nn, i, k
  real(wp) :: gr, dtn, dzn, dz, coeff1, coeff2
 
   nn = size(zn)
   gr = phc%gravity/phc%rgas

   do i=1,size(z)
     if(z(i).le.zn(1)) then
       p(i) = pn(1)
     elseif(z(i).gt.zn(nn)) then ! isothermal profile
       p(i) = pn(nn) * exp( -gr/tn(nn)*(z(i)-zn(nn)) )
     else
       do k=2,size(zn)
         if( z(i).le.zn(k) ) then
           dz  =  z(i)-zn(k-1)
           dzn = zn(k)-zn(k-1)
           dtn = tn(k)-tn(k-1)
           if(abs(dtn).gt.dt_toll) then
             coeff1 = tn(k)/tn(k-1)
             coeff2 = zn(k)-coeff1*zn(k-1)
             p(i) = pn(k-1) *                                     &
               ( ((coeff1-1.0_wp)*z(i)+coeff2)/dzn )**(-gr*dzn/dtn)
           else
             coeff1 = gr/tn(k-1)
             coeff2 = 0.5_wp*coeff1/tn(k-1)
             p(i) = pn(k-1) * exp( -coeff1*dz )              &
               * ( 1.0_wp + coeff2*dz**2/(zn(k)-zn(k-1))*dtn )
           endif
           exit
         endif
       enddo
     endif
   enddo
 
 end function pint

!-----------------------------------------------------------------------

end module mod_boulder_windstorm_test

